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Low-frequency, wide field-of-view (FoV) radio telescopes such as the Murchison 
Widefield Array (MWA) enable the ionosphere to be sampled at high spatial completeness. 

We present the results of the first power spectrum analysis of ionospheric fluctuations in 
MWA data, where we examined the position offsets of radio sources appearing in two 
datasets. The refractive shifts in the positions of celestial sources are proportional to 
spatial gradients in the electron column density transverse to the line of sight. These can 
be used to probe plasma structures and waves in the ionosphere. The regional (10-100 
km) scales probed by the MWA, determined by the size of its FoV and the spatial density 
of radio sources (typically thousands in a single FoV), complement the global (100-1000 
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km) scales of GPS studies and local (0.01-1 km) scales of radar scattering measurements. 
Our data exhibit a range of complex structures and waves. Some fluctuations have the 
characteristics of travelling ionospheric disturbances (TIDs), while others take the form of 
narrow, slowly-drifting bands aligned along the Earth’s magnetic field. 


1. Introduction 

1.1. Physical Background 

The ionosphere is the ionized component of the 
Earth’s atmosphere. It lies between ~50-1000km 
altitude, transitioning smoothly into the plasmas- 
phere, which is a torus of fully-ionized plasma ex¬ 
tending up to several Earth radii [Singh et at, 2011]. 
At night, the electron density peaks near 300-400 km 
[Luhmann, 1995; Solomon, 2010]. The ionosphere 
and plasmasphere are embedded with complex den¬ 
sity structures that can cause the refraction and 
diffraction of radio waves [Darrouzet et al., 2009; 
Schunk and Nagy, 2009]. These have been stud¬ 
ied for many decades with low-frequency radio inter¬ 
ferometer arrays [Hewish, 1951; Bolton et al., 1953; 
Mercier, 1986; Jacobson and Erickson, 1992a; Jacob¬ 
son et al., 1996; Helmboldt et al, 2012a; Helmboldt, 
2014]. 

Coupling between the various layers of the 
Earth’s neutral atmosphere and ionosphere oc¬ 
curs through electric fields, ion-neutral collisions, 
shear flows, and upward-propagating acoustic- 
gravity waves (AGWs). These mechanisms, as well 
as the effects of thunderstorms, geomagnetic sub¬ 
storms, solar X-ray flares and large seismic events 
can produce irregularities in the ionosphere [Rottger, 
1978; Luhmann, 1995; Lastovicka, 2006; Liu et al., 
2011]. Many travelling ionospheric disturbances 
(TIDs) are AGWs, which typically have wavelengths 
of 100-1000 km and periods of minutes to hours 
[Hocke and Schlegel, 1996]. These have been ob¬ 
served to produce quasi-periodic shifts in the posi¬ 
tions of radio sources as the waves propagate across 
the field of view [Bougeret, 1981]. 

Field-aligned density ducts are cylindrical en¬ 
hancements/depletions in the plasmasphere with 
widths of 10-100km [Angerami, 1970; Ohta et al., 
1996]. These can extend down to ionospheric alti¬ 
tudes [Bernhardt and Park, 1977] and may form by 
the interchange of magnetic flux tubes under the in- 
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fluence of electric fields, e.g. those from thunder¬ 
storms [Park, 1971; Walker, 1978]. They have been 
reported to produce rapid variations in interferom¬ 
eter visibility phases as they drift relative to the 
background sources [Jacobson and Erickson, 1993; 
Hoogeveen and Jacobson, 1997a; Helmboldt and In¬ 
terna, 2012]. 


1.2. Technical Background 


The apparent angular offset in the positions of ce¬ 
lestial radio sources (i.e. refraction) from their true 
coordinates results from the distortion of radio-wave 
phase fronts due to differences in the electron col¬ 
umn density (total electron content, TEG) between 
separate lines of sight to receivers [Thompson et al., 
2001; Schunk and Nagy, 2009]. The angular offset is 
given by 
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where e and rUe are the electron charge and mass, 
Co is the vacuum permittivity, n is the observing fre¬ 
quency, and VuTEG denotes the spatial gradient of 
the TEG transverse to the line of sight. The iy~^ de¬ 
pendence implies that offsets are larger at lower fre¬ 
quencies, and the negative sign indicates that sources 
refract towards lower TEG. The refractive wander¬ 
ing of a source over an integration due to temporal 
variations in the TEG can cause smearing and ap¬ 
parent decreases in peak flux density [Kassim et al., 
2007]. In addition, angular broadening and scintil¬ 
lation, which relate to the second spatial derivative 
of the TEG, occur when the angular sizes of density 
inhomogeneities are smaller than or comparable to 
the scattering angle [Lee and Jokipii, 1975a, b; Pick¬ 
ett, 1977; Meyer-Vernet, 1980; Meyer-Vernet et al., 
1981]. 

Methods for probing the ionosphere with radio 
waves include (1) ionospheric sounding, where the 
timing of broadband pulse reflections yields the ver¬ 
tical TEG profile [Munro, 1950; Reinisch, 1986; Lei 
et al., 2005]; (2) radar scattering, where powerful 
narrowband pulses reveal small-scale TEG irregular¬ 
ities [Gordon, 1958; Larsen et al, 2007]; (3) Global 
Positioning System (GPS) satellite measurements, 
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where dual-frequency delays yield the absolute TEC 
along given lines of sight [Mannucci et al., 1998; 
Lee et al., 2008]; (4) interferometer arrays, where 
TEC gradients are inferred from visibility phase fluc¬ 
tuations [Hewish, 1951; Hamaker, 1978]; and (5) 
satellite-borne beacon receivers, which measure the 
TEC based on differential delays of signals from the 
ground and other satellites [Bernhardt and Siefring, 
2006; Bernhardt et al, 2006]. Radar scattering can 
detect structures 1-10 m in size, over a region the 
width of the illuminating beam (several km). For 
GPS observations, the coverage and resolution de¬ 
pend on the distribution of ionospheric pierce points, 
which are the points where the lines of sight to 
ground stations intersect the ionosphere. Global 
GPS TEC maps with ~100km resolution are pub¬ 
licly available, e.g. from the Madrigal database 
[Rideout and Coster, 2006]. Instruments aboard 
satellites in low-Earth orbit receiving signals from 
ground stations, such as the Scintillation and TEC 
Receiver in Space (CITRIS) instrument, have en¬ 
abled the identification of scintillation-causing irreg¬ 
ularities on scale sizes between 100 m and 1 km, and 
TEC measurements at a spatial resolution of 7.5 km 
along the orbit path [Siefring et al, 2011]. However, 
these do not currently provide as global a coverage 
as GPS. 

Many radio interferometer studies of the iono¬ 
sphere have involved tracking a single bright source 
and measuring the phase fluctuations on each base¬ 
line (where a baseline refers to a pair of interferomet¬ 
ric receivers) as a function of time [Hamaker, 1978; 
Jaeobson and Erickson, 1992b; Coker et al, 2009; 
Helmholdt et al, 2012a, b]. The ionospheric phase (f) 
(oc A6) is proportional to the difference in TEC be¬ 
tween the lines of sight to the two receivers, so this 
can be used to probe temporal variations in the dif¬ 
ferential TEC on spatial scales corresponding to the 
projection of the array onto the ionosphere. We will 
refer to this approach as the antenna-based method. 

In antenna-based studies, the spatial sampling 
pattern is given by the array configuration, result¬ 
ing in direction-dependent spatial resolution and a 
poor impulse response for power spectrum analysis 
of spatial features with sparse arrays [e.g. Helmboldt 
and Interna, 2014]. Widefield instruments that are 
able to observe many celestial radio sources simulta¬ 
neously can overcome these limitations, because of 
the large number of isotropically-distributed pierce 
points. We will refer to the approach of analyzing 


source position shifts rather than antenna phases 
as the field-based method. There have been sev¬ 
eral field-based studies of the ionosphere that used 
the Very Large Array (VLA) at 74 MHz [Kassim 
et al, 1993]. These include the work of Cohen and 
Rottgering [2009], who conducted a statistical anal¬ 
ysis of pairwise differential offsets (3-8 sources in a 
FoV), and Helmholdt and Interna [2012], who per¬ 
formed a spectral analysis of the position offsets of 
29 bright sources in a single FoV and detected groups 
of wavelike disturbances. A combination of antenna- 
and field-based techniques has been used to conduct 
a climatological study of ionospheric irregularities us¬ 
ing the VLA, characterizing structure on scales of 
several to hundreds of kilometers [Helmboldt et al, 
2012c]. 

The Murchison Widefield Array [MWA; Lonsdale 
et al, 2009; Bowman et al, 2013; Tingay et al, 2013] 
is a 128-element, low-frequency (80-300 MHz) dipole 
array located at the Murchison Radio-astronomy Ob¬ 
servatory (MRO) in Western Australia. Each re¬ 
ceiver element, called a “tile”, comprises 16 dual¬ 
polarization bowtie dipole antennas arranged in a 
4x4 configuration. The 128 tiles are spread over 
a 3 km-wide area in a pseudo-random but centrally 
condensed configuration, yielding a very dense in¬ 
stantaneous u, u-coverage which gives it an excellent 
snapshot imaging capability. It has a ~30°-wide 
FoV, enabling it to image a patch of the ionosphere 
about 150 km across (assuming an ionospheric alti¬ 
tude of 300 km). Over 1000 unresolved point sources 
are typically visible in a single 2-min integration, 
yielding pierce-point separations of 1-10 km. With 
a high imaging cadence of up to 0.5 s (the cadence 
at which the raw visibilities are stored, which sets 
the minimum integration time of each snapshot), the 
MWA enables us to study ionospheric structures and 
dynamics on regional scales at high spatiotemporal 
resolution. This complements the global scales and 
lower spatial resolutions of GPS measurements, and 
the highly localized probing of radar scattering ex¬ 
periments. Airglow imaging [Makela et al, 2006; 
Shiokawa et al, 2013] can achieve widefield cover¬ 
age similar to the MWA, but its reliance on opti¬ 
cal wavelengths restricts operations to cloudless and 
moonless nights. Radio telescopes suffer no such lim¬ 
itations; radio observations can even be conducted 
in the daytime. Furthermore, airglow measurements 
are restricted to lower altitudes, whereas radio in¬ 
terferometers can detect irregularities over large al¬ 
titude ranges up to high altitudes (in the topside re- 
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gion and plasmasphere) which are otherwise difficult 
to access from the ground [Jacobson et al., 1996; Loi 
et al., 2015, in press]. 

The advent of sensitive, low-frequency, widefield 
radio instruments such as the MWA, the Low Fre¬ 
quency Array [LOFAR; van Haarlem et al., 2013] and 
the Long Wavelength Array [LWA; Ellingson et al., 
2009] gives us an opportunity to study the iono¬ 
sphere in unprecedented spatial detail and breadth. 
Characterization of the ionosphere is relevant not 
only to ionospheric science but also to radio astron¬ 
omy, because fluctuations in position and flux den¬ 
sity can complicate the cross-matching and classi¬ 
fication of radio sources. The ionosphere can also 
cause time- and position-dependent distortions of the 
telescope point-spread function (PSF), demanding 
non-trivial calibration and deconvolution schemes 
[Mitchell et al, 2008; Interna et al., 2009; Bernardi 
et al., 2011]. Information about the expected spatial 
and temporal scales of these fluctuations is important 
for designing and optimizing algorithms that identify 
and remove ionospheric distortions. 

1.3. This Work 

We present the results of a power spectrum anal¬ 
ysis of ionospheric fluctuations observed using the 
MWA. We analyzed the position offsets of point 
sources appearing in two MWA datasets from the 
commensal Epoch of Reionization (EoR) Transients 
Survey at 183 MHz (G0005; Pis: P. J. Hancock, 

L. Feng & N. Kudryavtseva) and the MWA Tran¬ 
sients Survey (MWATS) at 154 MHz (GOOOl; PI: 

M. Bell), hereafter datasets A and B respectively, as 
a function of sky position and time. 

This paper is structured as follows. In §2, we de¬ 
scribe the data reduction process and our approach 
to computing the power spectra. In §3, we demon¬ 
strate this on MWA data from the two aforemen¬ 
tioned datasets. We interpret the observed features 
in §4. Finally, we conclude in §5. 

2. Method 

All MWA data mentioned in this work were ob¬ 
tained at night. This work uses radio images syn¬ 
thesized using the entire array (field-based method), 
rather than information from individual baselines 
(antenna-based method). The whole MWA can be 
considered as a single receiver, and the number of 
pierce points is equal to the number of celestial ra¬ 


dio sources visible at any given epoch (e.g. Figure 
1). The antenna-based method, which needs to ob¬ 
serve a single bright source that dominates the vis¬ 
ibilities, is unsuitable for the MWA because of its 
large FoV and thereby large number of bright sources 
in any given field. Since the position offset of a source 
in an image formed using all 128 tiles is related to 
the average ionospheric phase over all baselines, the 
field-based approach effectively smooths structure in 
the ionosphere by a window of diameter equal to the 
physical dimensions of the MWA. The fairly compact 
size of the array (longest baseline ~3 km) means that 
our angular resolution for ionospheric work is about 
0?6, comparable to the average angular separation 
of point sources at typical instrument sensitivities. 
This is about an order of magnitude lower than the 
angular resolution of the images, which is set by the 
size of the MWA synthesized beam (~2-3 arcmin). 

2.1. MWA Data Reduction 

The process of synthesizing an image from inter¬ 
ferometric visibilities involves a Fourier inversion of 
the spatial coherence function of the electric field, 
which is sampled discretely by the individual base¬ 
lines [Thompson et at, 2001]. Prior to the inversion 


UTC 2013-09-16 19:00:39 



Right Ascension (J2000) 


Figure 1. Example snapshot image taken by the MWA 
(dataset B, epoch 53), zoomed in to the central one-third 
so that individual sources are visible. Each of the ~ 10^ 
celestial radio sources in each snapshot image contributes 
a pierce point. The intensity is displayed on a linear 
greyscale. 
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step it is often necessary to excise (“flag”) unusable 
portions of the data, for example those contaminated 
by radio-frequency interference (RFI), and also to 
calibrate for varying instrumental gains, which might 
arise for example from thermally-induced changes in 
cable lengths. Several algorithms exist for perform¬ 
ing instrument gain calibration. These include phase 
referencing [Fomalont and Perley, 1999], which in¬ 
volves observing a well-understood source of emission 
(a calibrator) and solving for the antenna gains that 
minimize the overall differences between measured 
and expected values, and self-calibration [Pearson 
and Readhead, 1984; Cornwell and Fomalont, 1999], 
which exploits a set of closure relations that relate 
the gains of each closed loop of antennas in an over¬ 
determined fashion to solve for individual antenna 
gains. 

Fourier transforming the gridded and weighted vis¬ 
ibilities produces an image (the “dirty” image) that is 
the true sky brightness distribution convolved with 
the telescope PSF. The shape of the PSF is deter¬ 
mined by the sampling pattern of the spatial coher¬ 
ence function (the u, u-coverage) and in general has 
a complicated shape. Various deconvolution algo¬ 
rithms exist for recovering the true brightness distri¬ 
bution from a dirty image. Point-like (unresolved) 
sources, such as those considered for analysis here, 
can be deconvolved effectively using the CLEAN al¬ 
gorithm [Flogbom, 1974]. This involves the iterative 
subtraction of a scaled version of the PSF from lo¬ 
cal intensity maxima, until certain convergence cri¬ 
teria are satisfied, to obtain a set of CLEAN compo¬ 
nents representing the point sources. These compo¬ 
nents are reconvolved with a two-dimensional Gaus¬ 
sian whose width is representative of the angular res¬ 
olution of the telescope and added back to the resid¬ 
uals of the deconvolution step, yielding an idealized 
approximation to the sky brightness distribution (the 
“restored” image). 

The two datasets discussed in this work were orig¬ 
inally obtained by the MWA for different scientific 
purposes (neither ionospheric), under different ob¬ 
serving strategies and reduced by separate science 
teams using different pipelines. We have appropri¬ 
ated these data for the current work; these differ¬ 
ences are not intentional and are beyond the scope 
of control of our study. 

Dataset A was taken on 2014-02-05 at 183 MHz 
with 30.72 MHz bandwidth, in “drift-and-shift” 
mode, in which the MWA analog beamformers point 


towards a new azimuth and elevation approximately 
every half hour, and a consistent patch of sky then 
drifts through the FoV. This minimizes the num¬ 
ber of pointings used, while still tracking the same 
patch of sky. Two-minute snapshots were continu¬ 
ously taken of the same patch of sky over the observ¬ 
ing duration, which was centred near local midnight 
(see Table 1 for details). 

Known failed antennas were flagged, comprising 
12% of the baselines. The data were flagged for 
RFI using the automated flagging routine AOFLAG- 
GER [Offringa et ai, 2012], removing 1% of the vis¬ 
ibilities. They were time-averaged from 0.5 s to 4 s 
and were not smoothed from their native 40 kHz 
frequency resolution. The data were calibrated 
from an observation of 3C444 (RA=22*'14™25.752®; 
Dec=—17°01'36'.'29), assuming a point-source model. 
This has a flux density of 67 Jy and a spectral index 
of—1.0 at 183MHz [Hurley-Walker et al., 2014], suf¬ 
ficient to dominate the visibilities and produce a set 
of per-channel, per-polarization, per-antenna com¬ 
plex gains. 

After applying these gains, the data were imaged 
using multi-frequency synthesis across all 30.72 MHz 
of bandwidth using wsclean [Offringa et ai, 2014], 
down to the first negative CLEAN component. This 
initial model was Fourier-transformed to produce a 
model visibility set. These were then applied in 
one round of self-calibration to obtain a single so¬ 
lution for the amplitude and phase of each antenna 
over the whole observing interval. This correction is 
position- and time-independent: it should have no 
effect on the relative positions of the field sources, 
and should only act to reduce instrumental errors. 
Note that unlike for conventional instruments with 
much narrower fields of view, solely antenna-based 
corrections are ineffective at removing ionospheric 
phases from MWA data, because its wide FoV im¬ 
plies that the ionosphere cannot be treated as lo¬ 
cally uniform within each FoV (and therefore that 
antenna and ionospheric phases are non-degenerate). 
The compact size of the MWA ensures that the ef¬ 
fective antenna-based ionospheric phasor, which is 
a visibility-weighted average of celestial source iono¬ 
spheric phasors over the FoV, is almost identical for 
each antenna thus allowing instrumental errors to 
dominate antenna-based phase differences. 

The data were then imaged in four 7.68 MHz sub¬ 
bands. The pixel resolution in all imaging steps was 
0'.75 and imaging included the entire FoV down to the 
first null of the primary beam (full width at first null 
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~50°). This yielded images of size 5000 x 5000 pixels. 
The images were produced in instrumental Stokes 
(the E-W dipole response, N-S dipole response, and 
the cross-terms). Using a model of the MWA pri¬ 
mary beam [Sutinjo et ai, 2014], these instrumental 
polarizations were then combined to produce Stokes I 
images for onward analysis. 

Dataset B was obtained on 2013-09-16 at 154 MHz 
over a 30.72 MHz bandwidth using drift-scan obser¬ 
vations at a declination of —55?0 on the meridian. A 
two-minute snapshot was obtained every six minutes 
at this declination for an entire night of observations 
(see Table 1 for details). These data were flagged 
for RFI (using AOFLAGGEr) and known problematic 
tiles. They were then averaged by a factor of two in 
time and four in frequency to yield a time resolution 
of 1 s and a frequency resolution of 160 kHz. 

Frequency-dependent, time-independent calibra¬ 
tion solutions as a function of tile and polarization 
were derived from an observation of PKS 2356—61 
using a source model extrapolated from 843 MHz 
[Mauch et al., 2003], and applied to all target ob¬ 
servations. This has a peak flux density of 30.6 Jy 
(integrated flux density of 114.6 Jy) at 154 MHz and 
a spectral index of —0.85. No self-calibration was 
performed on this dataset. The full-bandwidth data 
(30.72 MHz) were then imaged and deconvolved us¬ 
ing the WSCLEAN algorithm [Offringa et al, 2014] 
to yield images of size 3072 x 3072 pixels, the pixel 
size being 0'.75. The data were CLEANed up to a 
maximum of 20,000 iterations. The full 30.72 MHz 
bandwidth was imaged. For each observation, images 
for two instrumental polarizations (the N-S and E- 
W dipole responses) were corrected for the primary 
beam and then combined [for details, see Offringa 
et al., 2014] to produce a Stokes I image. 

2.2. Source Extraction and Selection 

We chose dataset A for discussion in this work be¬ 
cause the sources within it exhibited a moderately 
large root-mean-square (RMS) amplitude of position 
shifts, this being 19 arcsec on average and high com¬ 
pared to other datasets (characteristic RMS ampli¬ 
tude of 10-20arcsec at similar frequencies). This is 
smaller than the pixel size (45 arcsec) and the syn¬ 
thesized beam FWHM (119arcsec at 183MHz). We 
used the source finder Aegean [Haneock et al, 2012] 
to extract sources from the images, obtaining 1300- 
2100 sources per snapshot for a clip level of 5cr. Can¬ 
didate sources were identified by Aegean as local 


regions of negative spatial curvature whose maxima 
lay above a given noise threshold. Their parame¬ 
ters were then computed by fitting two-dimensional 
Gaussians to the flux density distribution at those 
locations. An estimate for the position error associ¬ 
ated with fitting a Gaussian to an unresolved source 
[Condon, 1997] is 

, ( 2 ) 

SNRy8hr2 

where 9b is the FWHM of the synthesized beam, as¬ 
sumed to be circular, and SNR refers to the signal-to- 
noise ratio, which we take to be the ratio of the fitted 
peak flux density to the local RMS noise. Note that 
positions can be measured to significantly greater 
precision than the instrument resolution, and that 
this precision increases with the SNR of the source 
(e.g. for a source at the 5 -ct detection threshold, 
its position can be measured correct to ~ 10 arcsec, 
an order of magnitude smaller than the synthe¬ 
sized beamwidth). We excluded candidate sources 
with negative fitted peak flux densities, and those 
with large position errors (>1 arcmin) as quoted by 
Aegean. When Aegean encounters an island of 
pixels for which there are fewer pixels than fit pa¬ 
rameters (6 parameters per component), the source 
parameters are estimated rather than fitted. We ex¬ 
cluded such sources from our analysis. 

Gross-matching the remainder with the NRAO 
VLA Sky Survey (NVSS) catalog [Condon et al., 
1998] using a cross-matching radius of 1.2 arcmin re¬ 
sulted in the identification of 4713 distinct sources 
and 55070 occurrences of any source in any epoch, 
implying that each source was observed in around 12 
epochs on average. This dataset contains 38 epochs 
in total. Cross-matching was performed by identi¬ 
fying the nearest NVSS source within the threshold 
radius of a given Aegean source and associating it to 
the latter. Sources extracted by Aegean for which 
there was no NVSS source within the search radius 
were discarded. Around 7% of sources were excluded 
in this manner. The purpose of cross-matching with 
an external catalog was to reduce the fraction of spu¬ 
rious sources, which may arise from noise spikes or 
residual sidelobes. 

We chose dataset B to contrast dataset A in 
that the average amplitude of position offsets was 
low (9 arcsec) compared to what is typically mea¬ 
sured in other datasets, about half the amplitude 
seen in dataset A. Source extraction returned 1600- 
2800 sources per snapshot for a clip level of 5 ct. 
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We subjected the sources to the same filtering and 
cross-matching criteria as for dataset A, except that 
we used the Sydney University Molonglo Sky Sur¬ 
vey (SUMSS) catalog [Mauch et al., 2003] for cross¬ 
matching to suit the relatively southern declination. 
A total of 19979 distinct sources were identified, and 
there were 179450 occurrences over all epochs. Each 
source was thus observed in around 9 epochs on av¬ 
erage. There are 77 epochs in total for this dataset. 
Sources for which there was no SUMSS match were 
discarded, this excluding ~ 1 % of sources. 

Differences between the two datasets in the num¬ 
ber of sources extracted per snapshot may be at¬ 
tributed to differing primary beam sizes (larger at 
lower frequencies), and differences in imaging and 
calibration approaches yielding different sensitivi¬ 
ties. Differences in the cross-matching efficiency be¬ 
tween the two datasets may be due to differing com¬ 
pleteness levels of the reference catalogs, or perhaps 
higher levels of ionospheric activity in dataset A 
causing sidelobe degradation and an increase in spu¬ 
rious source identifications compared to dataset B. 

2.3. Analysis of Position Offsets 

The power spectrum analysis technique described 
here is broadly applicable to any scalar quantity mea¬ 
surable over the pierce points. In this work we ana¬ 
lyzed just position offset information, where we con¬ 
sidered scalar projections of the position offset vector 
field onto different directions on the sky (we obtained 
these by taking the inner product of each vector with 
a unit vector pointing in a chosen direction). We 
computed the position offset vectors as the displace¬ 
ments of sources from their time-averaged celestial 
coordinates, which we took to be the reference points. 
The technique can also be applied to source flux den¬ 
sities, polarization vectors or rotation measures, but 
we defer these to future work. 

The purpose of using the time-averaged position of 
each source as a reference was to remove calibration- 
related position errors. We checked for these in 
separate analysis by cross-matching the sources in 
each dataset with the International Celestial Refer¬ 
ence Frame Second Realization (ICRF2) catalogue 
[Fey et ai, 2009]. Many datasets displayed collec¬ 
tive shifts of all sources from their ICRF2 positions, 
in some cases by up to one arcminute, where the 
direction of the shift depended on the calibrator. 
Since the calibration solutions have no time depen¬ 
dence, measuring the offsets with respect to the time- 


averaged position removes systematic position offsets 
arising from calibration errors, which are likely to 
produce spurious peaks in the power spectra. How¬ 
ever, the disadvantage of this operation is that iono¬ 
spheric fluctuations whose drift speeds happen to 
match those of the sources are subtracted away from 
the data. The use of a time average also alters the 
temporal spectral response, since this is akin to a 
detrending/high-pass filtering operation with an ef¬ 
fective window size given by the duration over which 
a source appears (~ 12 epochs = 24 min on average 
for dataset A, and ~9 epochs = 54 min for dataset 
B). The temporal frequency resolutions for the two 
datasets therefore differ from one another (roughly 
twice as high in dataset B than A), and are sys¬ 
tematically lower than given by the reciprocals of 
the total observing durations. However, since many 
sources are detected for longer durations than aver¬ 
age, the filtering effect is expected only to be partial 
and some sensitivity to lower temporal frequencies 
would be retained. 

We specified pierce-point locations in terms of a 
two-dimensional Cartesian coordinate system whose 
x-axis points west and y-axis points north. These 
were obtained via an equiareal projection of the ce¬ 
lestial coordinates of each radio source. The na¬ 
tive spherical coordinates for the projection [cf. Cal- 
abretta and Greisen, 2002], whose longitude and lat¬ 
itude we refer to as j5 and 7 respectively, are de¬ 
fined such that the fiducial point (/3o)7o) = ( 0 , 0 ) 
corresponds to the zenith and the South Celestial 
Pole lies at (/ 3 p, 7 p) = (90° — A,0°). Here A denotes 
the latitude of the observing site, measured positive 
southwards; for the MWA, A = 26?7. We have used 
the symbols (/3, 7 ) rather than the usual {(f), 9) of 
Calahretta and Greisen [2002] to avoid clashes with 
other variables defined in this paper. Geometrically, 
this corresponds to a 90° rotation of the Az/El coor¬ 
dinate system so that the coordinate equator lies on 
the meridian. 

The coordinates {x, y) were computed using the 
pseudocylindrical Sanson-Flamsteed projection 

0 ~ ((/3i-/3 )cos7-/3i) ’ 

where (/3, 7 ) = (/3i, 0) is the center of projection; here 
we chose /3i to be the median value of /3 over all pierce 
points in each dataset. With this choice of coordi¬ 
nates the a;-axis points west, the y-axis points north, 
and the origin [x, y) = (0, 0) lies at the zenith. The 
Sanson-Flamsteed projection preserves spatial scales, 
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although it fails to preserve angles. Angular distor¬ 
tions are minimized by our choice of native coordi¬ 
nates, whose poles are shifted away from the merid¬ 
ian where most MWA observations are conducted. A 
desirable property of these coordinates is that they 
are fixed to the terrestrial rather than the celestial 
sky, meaning that the velocity of a waveform will be 
measured with respect to the ground rather than the 
pierce points (which drift at ~20ms“^ through the 
ionosphere at an altitude of 300 km). 

For a given point source whose position offset vec¬ 
tor at some epoch is (Ax, Ay), we defined its projec¬ 
tion onto the 0-direction by 

S 0 = Ax cos 0 -I- Ay sin 0 , (4) 


respect to each of its three axes. We normalized the 
power spectrum by the observing volume V (size of 
the Fourier datacube, with units of deg“^ s“^) and 
the number of pierce points |Q| to obtain the dimen¬ 
sional power spectral density (units of deg"* s). This 
is 


'PimniO) = 


I 


exp 


Nx Ny Nt 

EEE Myfc(0) X 

i—1 j—1 k—1 

{i- l){l- 1) {j - l){m- 1) 


—27ri 


Nr 




{k — l)(n — 1) 
'Nt 


( 7 ) 


where the angle 0 is measured anticlockwise from the 
positive x-axis (0 can be thought of as the polar an¬ 
gle in 2D cylindrical coordinates). The scalar sg is 
the quantity analyzed in all subsequent sections. 


2.4. Power Spectrum Computation 

Every occurrence of a known NVSS or SUMSS 
source at a given epoch in a dataset induces a corre¬ 
sponding pierce point defined in terms of one tempo¬ 
ral (t) and two spatial (x, y) coordinates. The set of 
pierce points Q, along with their associated sg values 
computed as per Equation (4) for a chosen value of 
0, were then used to populate a three-dimensional, 
rectilinear, evenly-spaced datacube M according to 


Myfc(0) 


sg{q) q = {x,y,t) GQ 
0 otherwise. 


( 5 ) 


Here i,j,k € 1 are the array indices of M along the 
X, y and t directions, respectively. These were ob¬ 
tained from the spacetime coordinates of the pierce 
point q = (x, y, t) by gridding according to 

i{q) = [(X - Xinin)/d] 

ji.q) = \(.y-y^in)ld\ 

k{q) = t , (6) 


where N^, Ny and Nt are the dimensions of M along 
the X, y and t axes, respectively, i,m,n G Z+ index 
the kx, ky and uj dimensions of the output datacube 
V, and we have denoted t to avoid clashes 

with the symbols used for the indices of M. Here 
u) is the temporal frequency, kx is the E-W spatial 
frequency, and ky is the N-S spatial frequency. 

The temporal grid resolution for each dataset in 
dimensional units is simply equal to the snapshot 
cadence, which is 2 min for dataset A and 6 min for 
dataset B. To choose a spatial grid resolution, we 
considered the angular size of the MWA projected 
onto the ionosphere, since this represents the size of 
the window over which differential baseline phases 
are averaged. Given a longest baseline of about 
3 km and assuming an altitude of 300 km, this cor¬ 
responds to a resolution of 36 arcmin and a Nyquist 
sampling limit of 18 arcmin. We chose a spatial in¬ 
crement of d = 9 arcmin to grid the pierce points, 
which corresponds to spatially oversampling by a fac¬ 
tor of two. The sizes of the resulting datacubes were 
Nx X Ny X Nt = 264 x 229 x 38 elements for dataset 
A and 262 x 262 x 77 elements for dataset B. 

3. Results 


where d is the chosen grid resolution, Xmin and ymin 
are the minimum x and y values over Q, and t is the 
epoch number. Note that it is impossible to ensure 
that the mapping in (6) is unique. Where clashes 
occurred, we chose to average the sg values together 
(this occurred at a low rate of 2-3%). 

The spatiotemporal power spectrum V for each 
dataset was then calculated by taking the squared 
modulus of the discrete Fourier transform of M with 


Figure 2 shows the the position offset vector field 
for two selected snapshots from dataset A. There is 
a remarkable level of spatial organization into what 
appear to be patches of coordinated motion, where 
source displacements alternate between east and west 
on a spatial scale of ~5° between reversals (~I0° spa¬ 
tial periodicity). The patches are elongated roughly 
NW-SE (geographic). In the movie of the full dataset 
(see online supplementary material), these patches 
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Figure 2. Spatial distribution of position offset vectors 
for two selected snapshots in dataset A, where the mid¬ 
dle of each arrow marks the time-averaged location of a 
point source. The x-axis points west, the j/-axis points 
north, and (0,0) marks the zenith. Arrows are colored red 
if they have a negative x-component (eastward-pointing), 
blue if their r-component is positive (westward-pointing). 
Arrow lengths are scaled to 100 times the actual dis¬ 
placement distance. A movie showing the time lapse of 
this vector field is included as supplementary material. 
The data are plotted using a Sanson-Flamsteed projec¬ 
tion scheme (see §2.3 for details). See Movie SI for an 
animation of the vector field for the remaining snapshots 
of this dataset. 

can be seen to drift slowly across the sky towards 
the NE. 

In the case of dataset B, which was chosen for 
its relatively low level of ionospheric activity as in¬ 
ferred from the small magnitudes of position offsets 


Figure 3. As for Figure 2, but for dataset B and with ar¬ 
row lengths scaled at 120 times the actual displacement 
distance. See Movie S2 for an animation of the vector 
field for the remaining snapshots of this dataset. 

compared to other MWA datasets, the position offset 
vector field displays isolated N-S elongated structures 
in the form of narrow vertical streaks (see Figure 3 
for selected snapshots, and the corresponding movie 
in the supplementary material). These are much less 
pronounced than the bands in dataset A, and unlike 
in dataset A, there is no obvious periodicity in these 
structures. 

The degree of anisotropy of a vector field can be 
quantified by examining how the average projected 
amplitude of the vectors varies as a function of 0 . 
The total power of fluctuations in each dataset is de- 
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fined to be 


m = 


V 


Nx Ny 


Nt 


N^NyNt 


(8) 


l—l m—l n—1 


(integral of P(0) over all spatial and temporal fre¬ 
quencies, for a given projection angle 9). The de¬ 
pendence of this on 9 is plotted in Figure 4. We see 
that the fluctuation power per pierce point is higher 
in dataset A than B, which is consistent with the 
visibly larger amplitude of motion in dataset A. 

Figure 4 shows that T{9) maximizes along the E- 
W direction {9 fv 0°) in dataset A, and along the 
NW-SE direction {9 Ri 65°) in dataset B. The act 
of projecting the vector field onto a given value of 
9 is analogous to passing light through a polarizing 
filter, and so we can define the degree of anisotropy 
(cf. linear polarization) to be (M — m)/{M + to), 
where M and to are the maximum and minimum 
values of T(9), respectively. If the ionosphere is the 
main source of systematic anisotropies in the position 
offset vector field, then the degree of anisotropy is a 
nominal lower bound on the contribution of the iono¬ 
sphere to the angular position offsets. The degree of 
anisotropy in dataset A is ^30%, and in dataset B it 
is -10%. 

The power minima of Figure 4 correspond to char¬ 
acteristic residual position offsets of about ISarcsec 
for dataset A and 7 arcsec for dataset B. These are 
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Figure 4. The total power T{9) normalized to the num¬ 
ber of measurement points in each dataset, as a function 
of the orientation of the axis of projection 9, measured 
north of west. 


upper bounds on the contribution of isotropic sources 
of error (e.g. fitting errors, for roughly isotropic u, v- 
coverage) to the measured scatter in angular posi¬ 
tion. The mean value of the position fitting error A^/> 
(see Equation 2) is around 3 arcsec for both datasets, 
implying that fitting errors can only account for a 
small fraction of the position scatter. Assuming no 
other time-dependent effects, this indicates that a 
significant amount of the observed position shifts at 
154 and 183 MHz is due to the ionosphere. 

Eigure 5 shows a linear dependence of position off¬ 
sets on A^, where A is the observing wavelength, for 
bright (> 30 ct) sources. Such a dependence would 
be expected if the shifts were largely due to refrac¬ 
tive effects (cf. Equation (1)). However, the small 
fractional bandwidth of our observations prevents us 
from ruling out other functional dependencies. Al¬ 
though the dependence on A^ is well described by 
a straight line, the best-fit intercepts are non-zero. 
This suggests that multiple factors with different A- 
dependences may be contributing to the position off¬ 
sets, but a full investigation into this is beyond the 
scope of this paper. 

In the following subsections we present the power 
spectra computed for the two datasets over a range 



J? (m^) 


Figure 5. The magnitude of the angular position off¬ 
set of sources in four frequency sub-bands, as a function 
of squared wavelength. Each data point was obtained by 
averaging the magnitude of offset over all sources appear¬ 
ing in that band, restricted to those brighter than SOcr. 
The lines are linear fits to the four data points. The error 
bars on each data point associated with random fitting 
errors are smaller than the marker size. 
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of 9 values. The power spectra for several represen¬ 
tative values of 9 are included in this paper. For 
other 9 values, we refer the reader to supporting ma¬ 
terial comprising movies showing the evolution of the 
power spectra as the projection axis rotates through 
a full cycle. Since a unidirectional vector field van¬ 
ishes when the axis of projection is orthogonal to it, 
by rotating the axis of projection we selectively fil¬ 
ter out different wave modes and enhance others. In 
this manner, we infer the presence of multiple groups 
of waves propagating in a variety of directions and 
with different phase speeds. We verified with simu¬ 
lated data (see §A1) that our power spectra return 
results consistent with expectations. 

3.1. Response Function 

A rectilinear datacube of quasi-uniformly dis¬ 
tributed points would be expected to produce sinc- 
like ringing in both space and time; however, the 
natural sampling pattern of the MWA is influenced 
by two additional factors. Firstly, the MWA primary 
beam response gradually decreases towards the edges 
of the FoV, leading to a smooth decline in the S/N 
ratio and thereby a decline in the density of sources. 
This tends to soften spatial ringing. Secondly, the di¬ 
agonal tracks with respect to (x, y, t) made by pierce 
points drifting across the sky are expected to produce 
diagonal ringing features with respect to { kx , ky , ui ) 
in the response function, extending orthogonally to 
the drift direction in (x,y,t). Since celestial sources 
drift westward (increasing x with increasing t), we 
expect the ringing to be inclined towards negative 
kx with increasing to. 

The response functions for datasets A and B af¬ 
ter applying a Gaussian taper in the temporal di¬ 
mension are shown in Figs 6 and 7. These each 
have a compact main lobe with minimal side lobes 
(note the logarithmic scale and d-function-like spike 
at the origin), which is desirable for spectral anal¬ 
ysis, but there is evidence for the diagonal ringing 
patterns anticipated above. These features are more 
prominent in dataset A than B, and the main rea¬ 
son for this is the difference in survey mode be¬ 
tween EoR and MWATS. Because the EoR survey 
tracks a fixed patch of celestial sky, the spacetime 
distribution of pierce points is concentrated about 
an axis that is skewed westward following the di¬ 
rection of source motion. Longer on-source times 
also imply that source tracks extend for greater dis¬ 
tances on average than for MWATS. The drift speed 


of sources can in fact be calculated from Figure 
6b as the gradient of the diagonal features: this is 
duj/dkx ~ 4 X 10“^degs“^, which agrees with the 
angular rotation speed of the Earth at the latitude of 
the MWA (360° cos 26?7 day-1 = 3.7x degs"!). 

3.2. Power Spectra for Dataset A 

The power spectra for 9 = 45° and 135° are shown 
in Figs 8 and 9, respectively. The dominant mode 
for 9 = 45° is a low spatial-frequency mode with 
a period of ^Ihr, where the power spectral density 
distribution is elongated along the direction 50° anti¬ 
clockwise from the positive kx axis. This corresponds 
to large-scale structure (angular scale of order the 


(a) 



-1 - 0.5 0 0.5 1 

(deg”'') 


Figure 6. The response function of dataset A, where 
points have been gridded to a spatial resolution of 9 ar- 
cmin, averaged over (a) oj and (b) ky . 
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FoV) with phase fronts aligned roughly NE-SW trav¬ 
elling towards the NW. We discuss this in §4.2. 

For 9 = 135° there are strong peaks at |fc| « 
0.1deg“^, corresponding to phase fronts spaced reg¬ 
ularly by ^ 10° and stretching along 9 ~ 50° 
(i.e. NW-SE, geographic). This matches the direc¬ 
tion of the elongation seen for the bands in Figure 
2 and their spacing as estimated by eye. The lower 
panel in Figure 9 reveals the temporal frequency to 
be about 0.4mHz (period of 40min). The corre¬ 
sponding phase velocity is about 4 x 10“^ degs“^ to¬ 
wards the NE. There is additional power spreading 
out to higher spatial frequencies, indicating the pres¬ 
ence of smaller-scale structure up to |A:| ~ 0.5deg“^ 
(angular scales of ^ 2°). We discuss these observa¬ 
tions further in §4.1. 



Figure 7. As for Figure 6, but for dataset B. 


3.3. Power Spectra for Dataset B 

The power spectra for 0 = 0° and 90° are shown 
in Figs 10 and 11, respectively. In Figure 10, the 
power spectral density distribution is concentrated 
along ky — 0 over a substantial range of kx, out to 
\kx\ = 1 deg“^. This indicates the presence of small- 
scale structure with angular scales of ^ 1° highly 
elongated along the N-S direction, but with no strong 
periodicity. The lower panel (w, fca;-plane) shows no 
clear concentration of power associated with higher 
kx modes, meaning that these small-scale structures 
do not drift with a well-defined phase velocity. We 
discuss these in §4.1. 


0 = 45 ° 



-1 - 0.5 0 0.5 1 

(deg”'') 




Figure 8. Power spectrum averaged over (a) oj and (b) 
ky for dataset A, for 6 = 45°. See Movie S3 for an an¬ 
imation of the power spectrum for this dataset as the 
projection axis rotates. 
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There is a peak at w = 0 and the lowest non-zero 
kx value probed by the MWA that can be seen in Fig¬ 
ure 11, indicating the presence of a stationary mode 
with a wavelength of order the size of the FoV. We 
discuss possible explanations for this feature in §4.3. 

For 0 = 135° there appears to be a peak at around 
oj = 0.6 mHz near |A:| = 0, corresponding to a long- 
wavelength fluctuation of angular width of order the 
FoV with a period of about 30 min. However, we 
are unable to determine the propagation direction 
because this is not resolved in the fc^-plane. We 
discuss this in §4.2. 

3.4. Error analysis 

We performed an error analysis by comparing the 
amplitudes of the features in the power spectra with 


the noise level expected from position fitting er¬ 
rors. We used a Monte Carlo approach where we 
sampled position offsets from Gaussian distributions 
with characteristic widths given by Equation (2), and 
then computed the power spectra. We took the sig¬ 
nificance of the data to be the number of standard 
deviations above the mean, where the standard devi¬ 
ations and means were computed for each spatiotem- 
poral frequency over 200 Monte Carlo iterations. All 
features in the power spectra are many standard de¬ 
viations above the mean, indicating that the mea¬ 
sured signals are significantly larger than expected of 
random fluctuations. Figure 12 shows an example of 
a significance plot for dataset A, where the quantity 
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Figure 9. As for Figure 8, but with 6 — 135°. 


Figure 10. As for Figure 8, but for dataset B with 
9 — 0°. See Movie S4 for an animation of the power 
spectrum for this dataset as the projection axis rotates. 
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plotted is the number of standard deviations above 
the mean. 

An alternative means of estimating the error is to 
compare the measured positions between images gen¬ 
erated from the two instrumental polarizations (the 
N-S and E-W dipole responses), since these are effec¬ 
tively two independent, simultaneous measurements 
of the same signal. The ionosphere should have iden¬ 
tical effects on these images, and so any discrepan¬ 
cies represent the contribution of random noise. We 
find that the two corresponding sets of position mea¬ 
surements agree to within the position fitting errors 
predicted by Equation (2). An example comparison 
between the two polarization responses for a bright 
radio source is shown in Figure 13, where the posi¬ 


tion offsets for the two independent measurements 
are seen to be consistent within error bars. 

Although it is reassuring that the signals are 
strong compared to the amplitude of noise expected 
from fitting errors, the most problematic sources of 
error are likely to be those that are coherent rather 
than incoherent. The power spectral density of inco¬ 
herent noise distributes itself fairly evenly over the 
various frequency modes, and so should not interfere 
with the identification of coherent signals, which are 
localized in Fourier space. However, spurious struc¬ 
ture or the distortion of features in the power spectra 
can arise from a variety of causes, including sidelobes 
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Figure 11. As for Figure 10, but with 6 = 90°. 



Figure 12. Significance plot (showing number of stan¬ 
dard deviations above the mean) for the power spectrum 
of dataset A, averaged over (a) oj and (b) ky, for 9 — 90°, 
on a logarithmic scale. This was generated using 200 
Monte Carlo iterations. 
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of the response functions (see Figs 6 and 7), spectral 
leakage, and imaging artefacts. 

4. Discussion 

Similar power spectrum analyses of point-source 
position offsets have been performed previously by 
Helmholdt and Interna [2012] and Helmboldt et al. 
[2012c] using the VLA. However, the spatial sam¬ 
pling of our data is 1-2 orders of magnitude more 
complete and our FoV several times larger, which has 
not only yielded a better-behaved response function 
for power spectrum analysis, but has also allowed us 
to reconstruct ionospheric waveforms directly from 
the images. The MWA is the first radio telescope 
with the ability to do this; in previous studies making 



Time (min) 



Time (min) 


Figure 13. Comparison of the angular offsets of a bright 
(35-cr) source between images recorded in the two inde¬ 
pendent instrumental polarizations, along the (a) RA and 
(b) Dec directions (these are orthogonal). Error bars in¬ 
dicate the precision associated with Gaussian fitting. The 
two independent measurements are consistent with one 
another to within the measurement precision, which is 
significantly smaller than the amplitude of the position 
fluctuations themselves. 


use of radio interferometers [Jacobson and Erickson, 
1993; Jacobson et al, 1996; Hoogeveen and Jacobson, 
1997a, b; Helmboldt, 2014], the spatial sampling of 
the ionosphere was too sparse to visualize the den¬ 
sity structures in any detail. 

We have identified several types of fluctuations in 
the MWA data occupying different regions of Fourier 
space. These include the spatially-resolved, slowly- 
drifting fluctuations visible in Figs 2 and 3, and 
large-scale, unresolved fluctuations near jfc] = 0. Of 
the large-scale fluctuations, we have identified those 
with well-defined periods of 30 min and 1 hr, and also 
those that appear to be stationary with respect to the 
FoV (w = 0). 

4.1. Slowly-Drifting, Small-Scale Fluctuations 

As described in §3, spatially-resolved fluctuations 
were detected in both datasets A and B, where the 
structures appear to be elongated greatly along a cer¬ 
tain direction. Comparison with the Earth’s mag¬ 
netic field parameters at the location of the MWA 
(magnetic inclination of —60°, magnetic declination 
of 0°) reveals that these are consistent with be¬ 
ing irregularities that are extended along the geo¬ 
magnetic field lines. Their drifts where measured 
(4 X 10“^degs“^ in dataset A) are several times 
slower than would be expected of TID propaga¬ 
tion speeds at thermospheric heights (>100ms“^, or 
>0.02degs“^ at 300km altitude). This can be ex¬ 
plained if they are in-situ stationary structures that 
are drifting along with the background plasma. Drift 
speeds of several tens of meters per second are ex¬ 
pected for neutral wind motions at thermospheric 
heights [Amayenc and Vasseur, 1972]. 

At altitudes above about 150 km, where the elec¬ 
trical conductivity along the Earth’s magnetic field 
lines is substantially higher than transverse to them, 
the field lines can be treated as equipotentials [Park 
and Dejnakarintra, 1973]. At these heights, diffu¬ 
sion along the field lines is much faster than across 
them, and density enhancements that may initially 
form along a localized portion of the field lines read¬ 
ily diffuse along them to form large-scale, duct-like 
structures. These may then persist for significant 
durations, of order the timescales governing cross- 
field diffusion, which can be hours to days [Thomson, 
1978]. 

Although we have not estimated the altitudes 
of the structures seen here, previous analyses of 
field-aligned density ducts in MWA data by Loi 










16 


LOI ET AL.: MWA IONOSPHERIC POWER SPECTRUM ANALYSIS 


et al. [2015, in press] have obtained values of ~400- 
1000 km. If the structures lie at around 600 km 
altitude, their physical widths would be between 
10-100 km, similar to values measured by satellites 
for field-aligned ducts [Angerami, 1970] and the 
magnetic eastward-directed plasmaspheric waves de¬ 
scribed by Helmboldt and Interna [2012], The physi¬ 
cal drift speeds in dataset A would be about 40 ms“^ 
at this altitude. If the local fractional density varia¬ 
tions are comparable between the two datasets, the 
greater prominence of the structures in dataset A 
than dataset B can be explained if they are at a lower 
altitude in the former, since background densities are 
larger at lower altitudes. The suggestion that the 
structures in dataset A are nearer the ground, if the 
physical widths are intrinsically similar between the 
two datasets, is consistent with the larger angular 
widths observed in dataset A. 

Solar activity was significantly higher during the 
interval of dataset A (I0.7-cm radio flux index FI0.7 
= 188 sfu, where 1 sfu = 10“^^ Wm“^ Hz“^) com¬ 
pared to dataset B (F10.7 = 93 sfu). However, simi¬ 
larly extensive structuring is also seen in MWA data 
during periods of low solar activity, and so this ap¬ 
pears not to be a decisive influence on the observed 
behaviour. There is no clear link to geomagnetic ac¬ 
tivity that may explain the differences between the 
two datasets: both were obtained under quiet con¬ 
ditions where the Kp global storm index (quantify¬ 
ing fluctuations of the Earth’s magnetic field) was 
around 1, the 24-hr maximum Kp index not ex¬ 
ceeding 2 on either day. A broad study performed 
by Hoogeveen and Jacobson [1997a] using the Los 
Alamos radio interferometer found that field-aligned 
ducts of similar description appear most prominently 
on nights just after geomagnetic storms. Although 
we have yet to examine MWA data taken under sig¬ 
nificantly disturbed geomagnetic conditions, we have 
noted that field-aligned ducts appear routinely un¬ 
der quiet to mild conditions. Global GPS TEG maps 
[e.g. Rideout and Coster, 2006] have insufficient TEG 
precision and spatial resolution to reveal structures 
on these scales, but detections may still be possible 
with single-receiver relative TEC measurements. A 
climatological study of ionospheric behaviour above 
the observatory making use of combined MWA and 
on-site GPS TEC measurements is the subject of on¬ 
going work. 

The steeply inclined nature inferred for the density 
structures (60° to the ground, as given by the mag¬ 


netic inclination at the MRO) is at odds with the 
constant-altitude, thin-screen approximations com¬ 
monly used to model the ionosphere and employed 
in some calibration strategies [Rino, 1982; Wernik 
et al, 2007; Interna et al, 2009]. So far, no one 
attempting screen-based calibration has accounted 
for the possibility that density irregularities that af¬ 
fect radio observations may be localized to mag¬ 
netic shells rather than constant-altitude layers. The 
Earth’s magnetic field is well described by a dipole, 
and so the field lines are generally inclined with re¬ 
spect to the ground: they are horizontal at the mag¬ 
netic equator, vertical at the magnetic poles, and at 
intermediate inclinations at other latitudes. 

Although the bulk of the TEC is indeed concen¬ 
trated in a relatively thin layer at constant altitude, 
the constant offset component of the TEC is invisible 
to a radio interferometer, since this adds a constant 
phase to all antennas which then cancels out. In¬ 
terferometers differ fundamentally from conventional 
techniques for probing the ionosphere, the latter of 
which are mostly sensitive to absolute density and 
therefore tend to be dominated by the signal from 
the layer of peak electron density (300-400 km alti¬ 
tude). Because an interferometer is sensitive to 
TEC rather than absolute TEC, it experiences no 
such bias: irregularities at a large range of altitudes 
can easily dominate the signal, including those well 
above the peak density layer. 

The physically compact size of the MWA simpli¬ 
fies the problem of ionospheric calibration by mini¬ 
mizing the differences between the phase screens seen 
by different antennas [cf. Lonsdale, 2005], thus reduc¬ 
ing the dimensionality of the calibration problem to 
simply a function of angular position (besides time). 
However, screen calibration models that may eventu¬ 
ally be implemented on larger arrays (e.g. LOFAR) 
may need to allow for inclined screen components, 
to capture this tendency for magnetic shell localiza¬ 
tion. We find that duct-like, field-aligned (and there¬ 
fore steeply inclined) structures are very common in 
MWA data, appearing in half or more of all nights 
examined. This implies that a large population of ir¬ 
regularities causing measurable distortions in MWA 
data are ones that cannot be captured by a constant- 
altitude screen model. 

4.2. Large-Scale, Short-Period Fluctuations 

The long-wavelength (|fc| ~ 0), 30min-period fluc¬ 
tuations detected in dataset B are likely to be asso- 
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dated with a medium-scale TID (MSTID). MSTIDs 
have wavelengths of several hundred kilometers (of 
order the FoV or larger, at 300 km altitude) and pe¬ 
riods of several tens of minutes [Hunsucker, 1982], 
so these give rise to fluctuations occupying low- 
spatial frequency, high-temporal frequency regions 
of Fourier space. These are the most common type 
of ionospheric fluctuation believed to occur at mid¬ 
latitudes. We detect their spatiotemporal signatures 
in many other MWA datasets. 

The large-scale, time-resolved fluctuation in 
dataset A with a measured period of about 1 hr could 
be a large-scale TID (LSTID). LSTIDs have wave¬ 
lengths of ^1000 km and periods of an hour or more 
[Georges, 1968]. However, noting that its inferred 
propagation direction is parallel to the field lines, 
this may not in fact be a fluctuation confined to the 
thermosphere (as is the case for TIDs) but instead 
represent the wave-like transport of material up or 
down inclined magnetic flux tubes. While LSTIDs 
are expected to produce periodic, quasi-sinusoidal 
oscillations, the relatively short duration of our ob¬ 
servation for dataset A does not allow us to distin¬ 
guish between this and a one-off wave of deposition. 
Flows between the plasmasphere and ionosphere are 
expected to occur at night, when decay of the iono¬ 
sphere via recombination is combatted by downflows 
from the plasmaspheric reservoir [Park, 1970; Helm- 
boldt et al., 2012b]. 

There is a possibility that oscillating disturbances 
such as TIDs may be aliased at MWA survey ca¬ 
dences. The shortest TID wave period measurable 
at the cadence of the two MWA surveys used here is 
4 min for dataset A and 12 min for dataset B. Grav¬ 
ity waves, which comprise MSTIDs and LSTIDs, can 
only propagate at frequencies below the local Brunt- 
Vaisala (B-V) frequency, also known as the buoy¬ 
ancy frequency [Yeh and Liu, 1974]. The B-V fre¬ 
quency at ionospheric heights corresponds to a pe¬ 
riod of 10-12min [Hines, 1960], and so in the ab¬ 
sence of background winds these cadences should be 
adequate for characterizing MSTIDs and LSTIDs. 
However, background winds can Doppler-shift the 
waves to higher frequencies exceeding the Nyquist 
limit. The presence of winds with speeds as high 
as 100 ms“^ has been the proposed explanation for 
observations of gravity waves with periods as short 
as 5min [Davies et al, 1973]. Small-scale TIDs 
(SSTIDs), which are acoustic waves rather than grav¬ 
ity waves [Hunsucker, 1982], oscillate with frequen¬ 


cies above the B-V frequency and may be aliased at 
these cadences. However, raw interferometer visibili¬ 
ties are often stored at 0.5 s resolution, meaning that 
the data can be re-imaged under a much higher time 
resolution if necessary. 

4.3. Large-Scale, Stationary Fluctuations 

The peak at a; = 0 and the lowest non-zero spa¬ 
tial frequency in dataset B implies the existence of 
a fluctuation that is stationary over the FoV, with a 
length scale of order the FoV. This can be explained 
by way of either ionospheric or instrumental effects. 
Both are likely to contribute at some level, but we 
have not investigated into which may be the more 
dominant cause. 

A completely still, uniform ionosphere at a con¬ 
stant altitude can produce a systematic distortion 
in the angular position shifts of background sources 
because of the sec^ dependence of the TEC on 
the zenith angle C. The MWA FoV is about 30° 
wide, which implies that for a zenith pointing, the 
TEC along lines of sight near the edge of the FoV 
is about 4% larger than at zenith. For a back¬ 
ground TEC of order lOTECU (where ITECU = 
10^® electrons m“^) and an ionospheric altitude of 
300 km, the associated radial gradient is about 4 x 
10^*^ m“^, which would produce position shifts of or¬ 
der lOarcsec at 154 MHz. This is larger than typical 
position fitting errors and should be detectable via 
power spectrum analysis, since this signal would be 
coherent and persistent. Sources refract towards re¬ 
gions of lower TEC, and so the offset vectors associ¬ 
ated with this effect would be radially inward. 

Although this ionospheric distortion should be cir¬ 
cularly symmetric about the zenith, the signal ap¬ 
pearing in the power spectra for dataset B is con¬ 
centrated along E-W. This can be explained by the 
fact that the time-averaged position of each source 
is taken as the reference position. This causes sta¬ 
tionary features with respect to the sources to be 
subtracted out (recall the discussion in §2.3). Since 
celestial sources drift E-W, N-S components of the 
distortion pattern would be removed, leaving just an 
E-W signal. 

Instrumental effects that may give rise to large- 
scale distortions of source angular positions include 
tile position errors, and errors in correcting for the 
rc-terms during imaging [these describe the deviation 
of the array from coplanarity; Perley, 1999]. Unlike 
the secC distortion produced by the ionosphere, w- 
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correction errors can give rise to distortion patterns 
that are radially inward, outward, and/or with az¬ 
imuthal components. 

Applying a low-pass spatial filter to dataset B re¬ 
veals that the distortion pattern associated with the 
w = 0 peak is radially inward, potentially consis¬ 
tent with either source of error. However, in other 
MWATS datasets where w = 0 peaks appear, simi¬ 
lar checks sometimes reveal distortions that are ra¬ 
dially outward. This opposes the direction expected 
for ionospheric distortions, indicating that system¬ 
atic angular distortions may exist in MWA data that 
are not necessarily ionospheric in origin. 

4.4. Turbulence 

Kolmogorov turbulence [Kolmogorov, 1941] is ex¬ 
pected to produce a power spectrum that falls off 
as Note that although the 3D power spec¬ 
trum of TEC fluctuations goes as the quan¬ 

tity measured here is the TEC gradient. Taking a 
derivative multiplies the fluctuation spectrum by a 
factor of i\k\ and the power spectrum by a factor of 
leading to an overall dependence [Helm- 

boldt et al., 2012b]. The broad distribution of power 
spectral density along /cj, = 0 in the top panel of Fig¬ 
ure 10 may be evidence for a turbulent cascade, but 
this is not the only possibility. Narrow N-S density 
structures can also produce a power spectral density 
distribution that is broadly extended along ky = 0 
(cf. Fourier transform of a (5-function). That nar¬ 
row structures indeed appear in dataset B (Figure 3) 
suggests that these are the more likely explanation 
for the observed power along ky = 0, rather than an 
anisotropic turbulent spectrum, although we cannot 
rule out the existence of the latter. 

A further complication arises from the effects of 
spectral leakage and noise from incoherent measure¬ 
ment errors, both of which decrease the S/N ratio 
in the power spectra. These would flatten the power 
spectral density and preclude reliable measurement 
of the power-law index of a turbulent component, 
should there be one present (see §A2). A proper 
quantification of all major sources of error and de- 
convolution to remove sampling artefacts will be re¬ 
quired to measure the spatiotemporal properties of 
TEC fluctuations with greater accuracy. 

5. Conclusions &: Outlook 


We have performed the first power spectrum anal¬ 
ysis of ionospheric TEC variations using the MWA. 
Complex, coherent and highly organized fluctuations 
are often seen in the data. Power spectrum analysis 
is a useful technique for measuring the spatiotem¬ 
poral properties of these fluctuations. This comple¬ 
ments the superb ionospheric imaging capabilities of 
the MWA, which is the first radio telescope capable 
of imaging the ionosphere. Typical angular displace¬ 
ments at 154 and 183 MHz are 10-20 arcsec, sub-pixel 
most of the time but systematically larger than ex¬ 
pected from position fitting errors. The MWA is 
built on the planned site of the low-frequency por¬ 
tion of the future Square Kilometre Array [SKA-low; 
Dewdney et al, 2009]. It therefore has an impor¬ 
tant role to play in establishing the types and levels 
of ionospheric activity that will be experienced by 
SKA-low, and the analysis methods described here 
can be applied to this end. 

This work only investigated fluctuations in the an¬ 
gular position of point sources. However, flux density 
is another quantity of astronomical interest that may 
also be affected by the ionosphere. Ionospheric scin¬ 
tillation is expected to be induced by irregularities 
whose transverse spatial scales are several kilometers 
or less, which is below what we can directly mea¬ 
sure using the technique described here. The large- 
scale fluctuations we do detect should only produce 
negligible changes in the flux density, much smaller 
than the pixel-to-pixel RMS at typical instrument 
sensitivities. The current analysis therefore does not 
probe the TEC fluctuations that may give rise to sig¬ 
nificant amplitude variability. We intend to address 
this more quantitatively in a forthcoming paper. 

Our data show evidence for slow-moving, alternat¬ 
ing fronts of TEC depletion and enhancement which 
are likely to be field-aligned density ducts, and also 
longer-wavelength, faster-moving waves which may 
be TIDs. A stationary, radially-inward distortion 
pattern detected in dataset B may arise from a com¬ 
bination of ionospheric and instrumental effects. The 
presence of coherent, time-dependent fluctuations in 
both datasets, in particular dataset B which was 
chosen for its relatively low amplitude of source an¬ 
gular displacements, implies that ionospheric effects 
contribute significantly to source motions in MWA 
data at 154 and 183 MHz, and are present even in 
relatively undistorted datasets. This suggests that 
the successful implementation of algorithms target¬ 
ing ionospheric distortions will likely lead to signifi¬ 
cant improvements in astrometric accuracy for point- 
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source observations with the MWA and other low- 
frequency arrays. 

The wide instantaneous FoV and the high com¬ 
pleteness of pierce points enable the MWA to map 
ionospheric structure over a ~150 km region at a spa¬ 
tial resolution of ^3km (300km altitude), bridging 
the gap between global-scale GPS measurements and 
local-scale radar measurements. Contemporaneous 
GPS, radar, airglow and MWA observations would 
enable geophysical plasma phenomena to be stud¬ 
ied holistically on a wide range of length scales. Al¬ 
though airglow and MWA measurements probe sim¬ 
ilar spatial scales, the MWA is sensitive to a much 
larger range of vertical scales than airglow. In fu¬ 
ture, we intend to perform more extensive studies 
involving a larger number of MWA datasets, both to 
quantify the unwanted effects of the ionosphere on 
the primary data products for radio astronomy, and 
also to harness the valuable imaging capabilities of 
the MWA for the advancement of ionospheric science. 
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Appendix A: Simulated Data 

We performed experiments with simulated data to 
assess the severity of sampling-related artifacts re¬ 
lated to typical pierce-point patterns of the MWA. 
We tested two types of input: (1) a travelling plane 
wave, and (2) a spectrum of waves with amplitudes 
following a power law. To conduct the tests, we used 
the sampling pattern from dataset B, trimmed to di¬ 
mensions of 200 X 200 X 60 and with temporal gaps 
removed. 

Al. Experiment 1: Plane wave 

To test the system response we injected a single 
plane wave with a wavelength of 20 pixels and pe¬ 
riod of 4 epochs, travelling diagonally towards the 
bottom-right in Figure A2a. The resulting power 
spectrum is shown in Figs A2b and A2c, where we 
have collapsed the datacube down to two dimensions 
for visualization by averaging over either uj or ky. 

There are strong peaks at |fc| ~ 0.32 rad pix“^ 
and w Ri 1.6 rad epoch” implying a wavelength 
of ^20 pixels and a period of ^3.9 epochs, in close 
agreement with the input values. Note that the 
power spectra are inversion-symmetric and so a given 
wave mode produces a pair of peaks on opposite 
sides of the origin, corresponding to positive and 
negative frequencies of equal amplitude. The direc¬ 
tion of propagation can also be recovered from the 
power spectra: the line joining the peaks in Fig¬ 
ure A2b has dky/dkx = — 1, which matches the in¬ 
put propagation angle of —45°, modulo 180° (Fig¬ 
ure A2a). The wave is therefore either moving to¬ 
wards the top left or the bottom right. This can be 
narrowed down to a single direction by considering 
the line joining the peaks in Figure A2c, which has 
duj/dkx < 0 => dx/dt = —duj/dkx > 0, i.e. the 
phase velocity has a positive cc-component. Hence we 
conclude that motion is towards the bottom-right. 

A2. Experiment 2: Power law 

We tested a superposition of modes with ampli¬ 
tudes following a simple isotropic power law, where 
the power (squared modulus) of a certain mode with 
spatial and temporal frequencies of k and w, respec¬ 
tively, is given by 

V{w,k) rxk'^uj'^ . (Al) 

The input and output power spectra for m = 
— 11/3 and n = —\ are shown in Figure Al, inter¬ 
polated to a logarithmic grid. The flatness of the 
output spectrum compared to the input is likely to 


be a result of spectral leakage degrading the S/N ra¬ 
tio. However, the qualitative features of the input 
are recovered reasonably well. 

(a) 



-1 -0.5 0 0.5 

log^Q |k| (rad pix"'') 


Figure Al. (a) The input power spectrum for experi¬ 
ment 2, and (b) the output power spectrum, shown on 
logarithmic axes. The output spectrum has been boxcar- 
smoothed with a 3 X 3 pixel element with respect to the 
linear grid. Note the use of logarithmic axes, which are 
best for visualizing a power law. 
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Figure A2. (a) Snapshot of input plane wave for experi¬ 
ment 1, which has a wavelength of 20 pixels and a period 
of 4 epochs, and the output power spectrum averaged 
over (b) uj and (c) ky showing a peak at |fe| = 0.32 rad 
pix“^, UJ = 1.6 rad epoch“^. The diagonal ripples in (c) 
are features of the response function for dataset B, whose 
sampling distribution was borrowed for this experiment. 
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Table 1. Parameters for the two MWA datasets. 


Parameter 

Dataset A 

Dataset B 

Central frequency (MHz) 

183 

154 

Bandwidth (MHz) 

30.72 

30.72 

Polarization 

Stokes I 

Stokes I 

Calibrator 

3C444 

PKS 2356-61 

No. of snapshots 

38 

77 

Missing data 

None 

Snapshots 57 and 61 

Snapshot cadence (min) 

2 

6 

Integration time per snapshot (s) 

112 

112 

Span of observations (hr) 

1.3 

8 

UTC start 

2014-02-05 15:13:51 

2013-09-16 13:30:39 

UTC end 

2014-02-05 16:29:03 

2013-09-16 21:24:39 

AWST^tart 

2014-02-05 23:13:51 

2013-09-16 21:30:40 

AWST end 

2014-02-06 00:29:03 

2013-09-17 05:24:40 

LST start 

08:03:18 

21:00:00 

LST end 

09:18:43 

04:55:19 

Pointing center (J2000) 

RA, Dec = 150°,-24?8 

Dec = — 55?0 along meridian 

Pixel resolution (arcsec) 

45 

45 

Synthesized beam FWHM (arcsec) 

119 

130 

Visibility time resolution (s) 

4 

1 

Visibility frequency resolution (kHz) 

40 

160 


^ Australian Western Standard Time (UTC +8) 



